####This program was used in a BATCH mode
###R CMD BATCH "--args num1 num2 R P pp" Table4a.R
## Out-of-sample evaluation, rolling scheme
## It computes the simulations num1 to num2 (for parallelization)
## R is the size of the sample used for estimation
## P is the maximum size on which the moments are evaluated
## pp is the risk exposure of the VaR forecast
variable <- commandArgs(trailingOnly=TRUE)
ndeb <- as.numeric(variable[1])
nfin <- as.numeric(variable[2])
R <- as.numeric(variable[3])
p <- as.numeric(variable[4])
pp <- as.numeric(variable[5])


omega0=0.00001
alpha0=0.045
beta0=0.95
student0=3

library(fGarch)
library(MASS)

maxp=p

#Simulation d'un TGARCH
simu_garch_t<-function(nn){
  spec = garchSpec(model = list(omega=omega0, alpha = alpha0, beta = beta0,shape=student0), cond.dist = "std")
  #spec = garchSpec(model = list(omega=omega0, alpha = alpha0, beta = beta0))
  return(garchSim(spec, n = nn))
}





simu<-function(ns,r,p,pp){
  set.seed(10*ns)

  
  ##### Log-likelihood of a GARCH model ########
  garchLLH = function(x,parm) {
    omega = parm[1]; alpha = parm[2]; beta = parm[3]
    z =x; ltv = mean(z^2)
    # Use Filter Representation:
    e = omega + alpha * c(ltv, z[-length(x)]^2)
    h = filter(e, beta, "r", init =ltv )
    llh = sum(log(h)+z^2/h)/length(x)
    return(llh)
  }
  #Log likekihod with constraint  on the parameters
  vrais=function(x,b){
    bb=b;
    bb[1]=b[1]^2;
    bb[3]=exp(-b[2]^2-b[3]^2);
    bb[2]=exp(-b[2]^2)-bb[3];
    return(garchLLH(x,bb))
  }
    n=r+p
    ## We simulate a GARCH model of size r+maxp
    data=simu_garch_t(r+maxp)
    data=as.vector(data)
  
    ##############################
    ##### Rolling estimator ######
    rolling<-function(h){
      data_est=data[h:(r+h-1)] #size of estimation file is R
      
      #### QMLE ####
      #Initial conditions
      vol=var(data_est)
      adeb=0.045
      bdeb=0.95
      odeb=vol*(1-adeb-bdeb);
      tdeb=c(sqrt(odeb),sqrt(-log(adeb+bdeb)),sqrt(log((adeb+bdeb)/bdeb)))
      ## Optimisation procedure
      log_lik=function(param){vrais(data_est,param)}
      fit <-optim(tdeb,log_lik)
      estimpar=fit$par
      if (fit$convergence>0){estimpar=tdeb} # to avoid problems
      omega=estimpar[1]^2	
      beta=exp(-estimpar[2]^2-estimpar[3]^2);	
      alpha=exp(-estimpar[2]^2)-beta
      ###############
      
      #### We build the conditional var
      z =data[h:(r+h)]
      ltv = mean(z^2)
      e = omega + alpha * c(ltv, z[-length(z)]^2)
      h = filter(e, beta, "r", init =ltv )
      c3=filter(c(0,h[-length(z)]),beta,"r",init=0)
      c2=filter(c(0,z[-length(z)]^2),beta,"r",init=0)
      h=as.vector(h)
      #The derivative of log vol with respect to the parameters from t=1 to t=r+P
      der=cbind(1/h/(1-beta),c2/h,c3/h)
      der[1,]=c(0,0,0)
      
      ######Hit sequence and true score
      epst=as.vector(as.numeric((z[1:r]/sqrt(h[1:r]))))
      #Correction wrt the true score function
      dlog=der[1:r,]/2
      espder=apply(dlog,2,mean)
      Ahat=(qnorm(pp)*dnorm(qnorm(pp)))*espder
      Vhat=ginv(t(dlog)%*%dlog/r)/2
      corscore=t(Ahat)%*%Vhat%*%Ahat
      score=dlog*(epst^2-1)
      #Centered Hits
      It=1*(pnorm(epst)<pp)-pp
      #Calculation of E[dlog.*I_{t-h}]
      Itm1=c(NA,It[-length(It)])
      Itm2=c(NA,Itm1[-length(It)])
      Itm3=c(NA,Itm2[-length(It)])
      espder1=c(mean(Itm1*dlog[,1],na.rm=TRUE),mean(Itm1*dlog[,2],na.rm=TRUE),mean(Itm1*dlog[,3],na.rm=TRUE))
      Ahat1=(qnorm(pp)*dnorm(qnorm(pp)))*espder1
      corscore1=t(Ahat1)%*%Vhat%*%Ahat1
      espder2=c(mean(Itm2*dlog[,1],na.rm=TRUE),mean(Itm2*dlog[,2],na.rm=TRUE),mean(Itm2*dlog[,3],na.rm=TRUE))
      Ahat2=(qnorm(pp)*dnorm(qnorm(pp)))*espder2
      corscore2=t(Ahat2)%*%Vhat%*%Ahat2
      espder3=c(mean(Itm3*dlog[,1],na.rm=TRUE),mean(Itm3*dlog[,2],na.rm=TRUE),mean(Itm3*dlog[,3],na.rm=TRUE))
      Ahat3=(qnorm(pp)*dnorm(qnorm(pp)))*espder3
      corscore3=t(Ahat3)%*%Vhat%*%Ahat3
      
      ###VaR forecasts from R+1
      volf=h[r+1] # Variance forecast
      epstf=z[r+1]/sqrt(volf) # ratio of real return to vol forecast
      Itf=1*(pnorm(epstf)<pp)-pp
      
      ##### et in the paper
      etf=Itf+qnorm(pp)*dnorm(qnorm(pp))/2*(epstf^2-1)
      scoref=(der[r+1,]/2)*(epstf^2-1)
      ##### etstar in the paper
      etstarf=Itf+as.vector(t(Ahat)%*%Vhat%*%(scoref))
      
      ###Normalization to have unit variance
      etf=etf/sqrt((pp*(1-pp)-(qnorm(pp)*dnorm(qnorm(pp)))^2/2))
      etstf=etstarf/sqrt((pp*(1-pp)-corscore))
      if ((p/r)<1){poids=-(p/r)^2/3}else{poids=(2*r)/3/p-1}
      Itcor=Itf/sqrt(pp*(1-pp)+poids*corscore)
      
      return(cbind(Itf,etf,etstf,Itcor,Itf/(pp*(1-pp)*pp*(1-pp)+poids*corscore1)^{0.25},Itf/sqrt(sqrt(pp*(1-pp)*pp*(1-pp)+poids*corscore2)),Itf/sqrt(sqrt(pp*(1-pp)*pp*(1-pp)+poids*corscore3))))
    }
    matrol=t(sapply(seq(1,p,1),rolling))
    
    ##### Test rolling - unconditional ######
    testIt=p*(mean(matrol[,1])^2)/(pp*(1-pp))
    testet=p*(mean(matrol[,2])^2)
    testetst=p*(mean(matrol[,3])^2)
    testItcor=p*(mean(matrol[,4])^2)
    
    #### Test rolling - covariance #####
    statcov<-function(etfop){
      OP=length(etfop)
      etm1=c(NA,etfop[-length(etfop)])
      etm2=c(NA,etm1[-length(etfop)])
      etm3=c(NA,etm2[-length(etfop)])
      xij1=(OP-1)*mean(etm1*etfop,na.rm=TRUE)^2
      xij2=(OP-2)*mean(etm2*etfop,na.rm=TRUE)^2
      xij3=(OP-3)*mean(etm3*etfop,na.rm=TRUE)^2
      mj3=etfop*(etm1+2*etm2/3+etm3/3)
      v3=sum(seq(1,3,1)^2/9)
      xiw3=(OP-3)*mean(mj3,na.rm=TRUE)^2/v3
      return(c(xij1,xij2,xij3,xiw3))
    }
    #Joint moments e_t.e_{t-h}
    xicov1=statcov(matrol[,2])  
    #Joint moments etst_tetst_{t-h}
    xicov2=statcov(matrol[,3])
    #Joint moments ITI_{t-h} corrected
    ets=matrol[,5]
    etm1=c(NA,ets[-length(ets)])
    xim1=(p-1)*mean(etm1*ets,na.rm=TRUE)^2
    ets=matrol[,6]
    etm1=c(NA,ets[-length(ets)])
    etm2=c(NA,etm1[-length(etm1)])
    xim2=(p-2)*mean(etm2*ets,na.rm=TRUE)^2
    ets=matrol[,7]
    etm1=c(NA,ets[-length(ets)])
    etm2=c(NA,etm1[-length(etm1)])
    etm3=c(NA,etm2[-length(etm2)])
    xim3=(p-3)*mean(etm3*ets,na.rm=TRUE)^2
    xicov3=c(xim1,xim2,xim3)
    return(c(ns,r,p,pp,testIt,testet,testetst,testItcor,xicov1,xicov2,xicov3))
  }

### We launch the program and display the results
options(width=1000)
toto=(matrix(sapply(ndeb:nfin,simu,r=R,p=p,pp=pp),nrow=19))
noquote(formatC(t(toto),digits=3,format="f"))
